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A UNIFIED THEORY OF NON-IDEAL GAS LATTICE BOLTZMANN MODELS 


LI-SHI LUO* 


Abstract. A non-idcal gas lattice Boltzmann model is directly derived, in an a priori fashion, from the 
Enskog equation for dense gases. The model is rigorously obtained by a systematic procedure to discretize 
the Enskog equation (in the presence of an external force) in both phase space and time. The lattice 
Boltzmann model derived here is thermodynamically consistent and is free of the defects which exist in 
previous lattice Boltzmann models for non-idcal gases. The existing lattice Boltzmann models for non- ideal 
gases are analyzed and compared with the model derived here. 

Key words, lattice Boltzmann method, kinetic theory, Enskog equation, phase transitions, multi-phase 
fluids, fluid dynamics 
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1. Introduction. In recent years, there has been significant progress in the development of the lattice 
Boltzmann equation (LBE) method [1, 2, 3, 4, 5], a novel technique developed for modeling complex systems. 
One particular application of the lattice Boltzmann method which has attracted considerable attention is the 
modeling of inhomogeneous fluids, such as multi-phase or multi-component fluids [6, 7, 8, 9, 10]. These flows 
are important, but are difficult to simulate by conventional techniques of solving the Navicr- Stokes equations. 
The main difficulty conventional techniques face is the existence of interfaces in inhomogeneous flow. There 
is ample evidence that the lattice Boltzmann models based on mesoscopic theory are particularly suitable 
for these systems [6, 7, 8, 9, 10]. There are fundamental reasons for the success of the LBE models. Besides 
their broad applicability, the LBE models can also serve as new paradigms in nonequilibrium statistical 
mechanics, much like the Ising model in equilibrium statistical mechanics. Many hydrodynamic systems 
far from equilibrium are difficult to simulate by using the Boltzmann equation directly. The LBE method 
provides a novel and efficient means to simulate systems far from equilibrium. The LBE models do not start 
at the macroscopic level; instead, they start at the mesoscopic level at which one can freely use a “ potential 
to model interactions in the system. Macroscopic or hydrodynamic effects naturally emerge from mesoscopic 
dynamics, provided that the mesoscopic dynamics possess the correct and necessary conservation laws and 
associated symmetries. 

Historically, the lattice Boltzmann equation was first developed empirically [1, 2, 3] from its predecessor 

the lattice-gas automata [11, 12]. This empiricism influences even the most recent lattice Boltzmann 
models [6, 7, 8, 9, 10]. Empirical lattice Boltzmann models usually have some inherent artifacts which are 
not yet fully understood. One particular problem with multi-phase or multi-component lattice Boltzmann 
models is the thermodynamic inconsistency: the equilibrium state in these models cannot be described by 
thermodynamics [8, 9]. Although this issue has been raised previously [8, 9], no progress has been made in 
solving this problem, despite its paramount importance. 

It has been recently demonstrated [13, 14] that the lattice Boltzmann equation can be directly derived 
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from the continuous Boltzmann equation. The method of Refs. [13, 14] is a general procedure to construct 
the lattice Boltzmann models in a systematic and a priori fashion. Through this procedure we can better 
understand the approximation made in the lattice Boltzmann equation. In this paper, we apply the method 
of Refs. [13, 14] to analyze the lattice Boltzmann equation for multi- phase fluids with non-ideal gas equation 
of state. We derive the lattice Boltzmann equation from the Enskog equation for dense gas in the presence 
of an external force. We obtain a lattice Boltzmann equation for isothermal multi- phase fluids which has 
the required thermodynamic consistency. In addition, we compare our model with the existing ones. 

2. Enskog Equation and LBE Model for Non- Ideal Gases. It is well known that the original 
Boltzmann equation only describes rarefied gases; it does not describe dense gases or liquids. In the Boltz- 
mann gas limit (BGL), N — > oo, m — > 0, and r — > 0, where N, m, and r are the particle number, particle 
mass, and interaction range respectively, and Nm — > finite, Nr 2 — ► finite, and Nr 3 — » 0. Thus, in the 
BGL, the mean free path l ~ 1/Nr 2 remains constant, while the total interaction volume Nr 3 goes to zero. 
Therefore, in the strict thermodynamic sense, the Boltzmann equation only retains the thermodynamic prop- 
erties of a perfect gas there is no contribution to the transport of molecular properties from intcrparticlc 
forces, although collisions influenced by interparticlc interaction are considered. In order to properly describe 
non-ideal dense gases, the effect of finite particle size must be explicitly considered. It was Enskog who first 
extended the Boltzmann equation to dense gases by including the volume exclusion effect [15], which leads to 
a non-ideal gas equation of state. The Enskog equation [15, 16, 17] explicitly includes the radius of colliding 
particles, r 0 , in the collision integral: 

(2.1a) d t f + Z-Vf + a.Vtf = J, 

(2.1b) J = Ji dp x [g(x + r Q r)f(x. £')/(x + 2r 0 r, £\) - g(x - r 0 r)f(x, £)f(x - 2 r 0 r, £,)] , 

where / is the single particle (mass) distribution function, £ and a are particle velocity and acceleration, g is 
the radial distribution function, r is the unit vector in the direction from the center of the second particle of 
f( x > £ 1 ) t0 ^c center of the first particle of f(x , £) at the instant of contact during a collision, and pt x is the 
collisional space of the second particle of /(x, £ x ). If we expand the collision operator J in a Taylor series 
about x, use the BGK approximation [18, 17, 19], and assume the fluid to be isothermal and incompressible 
[20], we have: 

(2.2a) 0 t / + CV/ + a-V e / = -|[/- /“”] + J ' , 

(2.2b) J' = -/«•> bpg(t- u) V In (p 2 g ) , 

where A is the relaxation time and / (0) is the local Maxwell equilibrium distribution function given by 
(2.3) f i0) = p(2n0)- D / 2 cxp[-(Z - u) 2 /20] , 

where D is the dimension of the £ space; p , it, and 6 = kB r j'fm are the mass density, the macroscopic 
velocity, and the normalized temperature (per unit mass); and k& and T arc the Boltzmann constant and 
temperature. The additional collision term in Eqs. (2.2), J', describes the volume exclusion effect [20], where 
9 = g{bp ), and b is the second virial coefficient in the virial expansion of the equation of state. It is assumed 
that the acceleration a is due to an external potential U(x) (per unit mass): a = -VU. 

If the acceleration a is assumed to be a constant within the distance travel over a small time interval 
St, then a formal solution of Eq. (2.2) can be obtained by integrating along a characteristic line £ over the 
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time interval 6 t : 


(2.4) f(x + £6t + ia£ ( 2 , £ + a $t, t + 5 t ) = e Sl9/X f(x, i, t) 

+ l e ~s t g/x f ‘ e t'g/x f( o)( x + £ t > + I at ' 2 i £ + at', t + t') dt' 

^ JO 2 

+e~ St9/x [ 1 e t>9/x J'(x + it' + \at' 2 , i + at', t + t') dt ' . 

Jo ^ 

If wc assume that S t is small enough and both f m and / arc smooth enough in phase space, we can neglect 
the terms of order 0(5?) or smaller in the Taylor expansion of Eq. (2.4), and obtain [13, 14, 20]: 

( 2 . 5 ) f{x + is t , i, t+5 t )~ fix, i, t) = -2- [fix, i, t) - f m ix, i, t)} 

+ J'(®, £, t) S t - a-V^/(x, £, t) S t , 

where r = X/St is the dimensionless relaxation time. It is obvious that the accuracy of the above equation 
is only first-order in time [0(S t )). Consequently, the accuracy of the lattice Boltzmann models derived from 
the above equation is also first-order in time at best. 

For isothermal fluids, the equilibrium distribution function can be obtained by truncation of the Taylor 
expansion of f (0) up to second-order in u. 


where 

(2.7) L 0 (i) = (27 T0)" D/2 exp (-£ 2 /2 6) . 

The phase space discretization has to be done in such way that not only all the hydrodynamic moments, 
but also their fluxes arc preserved exactly. This is accomplished by using Gaussian quadrature to compute 
the moments [13, 14], 

Following the procedure described in Refs. [13, 14], we can obtain the LBE models in both 2D and 3D 
lattice space [13, 14]. We use the 2D nine-bit model as a concrete example here. In this case, wc have the 
following equilibrium distribution function [13, 14]: 


(2.8) 

f ( ^ = w a p l + ^r 

u) 9(e Q • «) 2 
2c 4 

3u 2 

2c 2 

where 

f 4/9, 

a = 0, 


(2.9) 

e 

p 

ii 

M 

ID 

a= 1, 2, 3, 4, 



{ 1/36, 

a = 5, 6, 7, 8, 



(0, 0), a = 0, 

(2.10) e a = (cos (f) a , sin<£ a )c, a — 1, 2, 3, 4, 

(cos<£ a , sin 0 a ) \/2c, a = 5, 6, 7, 8, 

</> a = (a- 1 ) 7 t /2 for a = 1 4, and <p a = (a - 5)7 t/ 2 + 7r/4 for a = 5 8, and c = S x /S t = \f30, and S x is the 
lattice constant. Note that 9 is a constant here. 


= 1 


(t-u ) 2 

2 0 2 
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The forcing term, a , is unknown but it can be written in terms of an expansion in £ as follows: 

( 2 - n ) a-W = M«[c (0) +ci 1) 6+c;^ + -”] • 

If the above expansion is truncated, the first few coefficient can be easily obtained by using the 

following moment constraints: 

(2.12a) Jd£aV(.f = 0, 

(212b) Jd££aV{f = —pa, 

( 2 - 12c ) jdUiij a-V € / = -piaiu, + a jUl ) . 

Therefore, up to the order of O(u) and 0(£ 2 ), we have 


(2-13) a-^f — -pio(£)0 1 [(£ — u)+0 

Note that in the above expansion, only terms up to first order in u have been retained, because there is a 
overall factor of St in the forcing term, as indicated in Eq. (2.5), and both St and u arc small parameters 
of the same order in the Chapman-Enskog analysis of the lattice Boltzmann equation [20, 21, 22]. There 
arc other methods to compute the forcing term [20]. It shoulc: be stressed that every term in the Enskog 
equation must be treated equally to maintain the same order o accuracy. Specifically, the expansion of the 
forcing term must be of second order in £ and of first order in u, in order to be consistent with the expansion 
of the equilibrium distribution function, given by Eq. (2.8). 

Following the same discretization procedure for the equilibrium distribution function, we obtain the 
forcing for the 9- bit model 


(2.14) 


F a 


-3 w a p 



u) -j- 3 


(f q-tt) 

C 4 



a . 


The above forcing term satisfies the discrete counterpart of Eqs. (2.12). If only the first two moment equations 
m Eqs. (2.12) arc satisfied, and the third constraint of Eq. (2.12c) is replaced by e Q(J e QJ F a = 0 in the 
discrete case, then, the forcing term reduces to F a = — 3w a pc 2 e a - a . This is the forcing term often used 
in the literature [21, 22]. 

The additional collision term J given by Eq. (2.2b) can be explicitly written in the discrete form: 


(2.15) 


J' a = -/^ eq) bpg(e a - u)-V n (p 2 g) . 


Including the discretized J\ the lattice Boltzmann equation obtained is: 

(2.16) f a (x + e a 6 t , t + 6 t )~ f a {x, t) = -1 [/ Q (x, t) - /^ eq) (x, t)j g 

~bpgf^ ] {x, t) (e a - u)-V(p 2 g)S t -F a S t . 
The Navicr-Stokes equations derived from the above LBE mode are [20]: 

(2.17a) d t p + V • (pu) = 0 , 

(2.17b) 3,u + uVu=--VP + i/V 2 u + a, 

P 
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where the viscosity 


(2 r-g)8l 


and the pressure (or the equation of state) is given by 


(2.19) 


P = pO(l Pbpg). 


Obviously, the above is a non-ideal gas equation of state. For ideal gases such that 6=0 and g = 1, P 
and v reduce to previous results for ideal gases. The dependence of the viscosity v on g can be removed by 
replacing g in the BGK collision term by 1. 

Given the equation of state, the Helmholtz free energy density can be obtained as: 


( 2 . 20 ) 


V>(p) = 


p j --dp = p9 In p+bjgt 


That is, with either P or 0 given, one can derive all the relevant thermodynamic quantities from the free 
energy function 0. With the free energy and the equation of state defined, the Maxwell construction [23] to 
determine the co-existence curve becomes physically meaningful and consistent. The phenomenon of liquid- 
gas phase transition can be simulated using this model by changing the value of b J g dp (by adjusting b or 
g) in the free energy density 0 relative to the temperature 0 as indicated by Eq. (2.20). 

3. Comparison with Existing Models. A comparison with the existing models [6, 7, 8, 9, 10] is now 
in order. In the Shan and Chen model [6, 7], an arbitrary potential U(x) = U(p(x)) is explicitly given, and 
the change of velocity u due to U(x ) is given by 


Su = —VU (x) rS t = a rS t 


By substituting u with u + Su into the equilibrium distribution function, we have 


(3.1) f ( ^ ) = w Q p 1 + 


+ Hr + 


3(e Q *ix) 9(e Q -u) 2 3v? 




— 3 w a p -(e a -u) + 3 


e a ■ xi) 


In the above result, the first part is the usual equilibrium distribution which has an ideal gas equation of 
state built in. The second part is supposed to account for interaction or non-ideal gas effects, which leads 
to the identical forcing term given by Eq. (2.14). By combining the forcing term with the pressure in the 
Navier-Stokcs equation, i.e., V pO + W = V(p0 + U), the equation of the state becomes P = [pO + U(p)}. 
Thus, the non-ideal gas effects arc effectively mimicked by the potential U. Of course, the physical concept 
of this approach is incorrect and the immediate shortcoming is that the heat flux, and hence the energy 
balance equation, is incorrect [20]. Furthermore, the third part in Eq. (3.1), which is proportional to 5} and 
nonlinear in a, is not consistent with what is obtained from Eq. (2.4). 

We should also discuss a recent revision of the Shan and Chen model [10] in which a forcing term 
proportional to (e a - u) * F5 t is derived with some crude approximations: the force F ~ -W - 
bpOgV ln(p 2 p), and V accounts for the attractive part in the interaction. This model produces a non-ideal 
gas equation of state, P = pG(l + bpg) + V, as expected. However, the derivation of this model closely 
follows the derivation of the previous model. Therefore, these two models share the same problems, such as 
incorrect heat transfer. 
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A comparison with the model proposed in [8, 9] is slightly more elaborate. Stressing the consistency of 
thermodynamics and being inspired by Cahn-Hilliard’s model ( !4], Swift et al. [8, 9] start with a free energy 
functional, 

(3.2) *= jdx^WVpf + Tfiip)] , 

where i}> is the bulk free energy density. The free energy functional in turn determines the diagonal term of 
the pressure tensor: 

(3-3) P = p^- - ^ = p- K pV 2 p - ^||Vp|| 2 , 

where p = pi’' - ip is the equation of state of the fluid. The full pressure tensor is given by 

(3-4) Pij — P 6 r j + Kdipdjp . 

With P tJ given, the equilibrium distribution function, / c ( , eq) , is constructed by not only satisfying the con- 
servation constraints, but also producing the above pressure tensor by forcing the following constraint: 

( 3 -5) £/ir°e a , <Ca j = p tJ . 

a 

It should be pointed out that in the context of Chapman-Enskog analysis, the presence of Vp related 
terms in 4' and P t l is not justified at all the density gradient does not appear in the first order Chapman- 
Enskog solution. Also, the model produces a number of unphysical effects. First, the term related to 
non-ideal gas effect misses a factor of (e a - u), and is therefore not Galilean invariant, as previously noticed 
[8, 9, 10], Second, the term related to e a e a , denoted as G tJ e a l e QiJ in Refs. [8, 9], is anisotropic, because 
G xx — -G yy . Third, the ratio between the number of the rest particles and the number of moving particles 
depends on the local density gradient. It can be shown that this ratio is directly related to temperature [25]. 
While the model is supposed to simulate an isothermal fluid, the temperature in this model may vary locally 
depending on the density gradient. Also, the model cannot lead to the correct energy balance equation. 
Furthermore, the pressure tensor P tJ does not appear in the Na vier-Stokes equation derived from the model 
[8, 9]. Therefore, the approach in deriving this model in Refs. 8, 9] is not only mathematically ad hoc and 
inconsistent, but also physically incorrect. 

It should also be noted that the Hamiltonian approach [6, 7] and the free energy approach [8, 9] arc 
indeed equivalent. Given the Hamiltonian of an interacting A'-particlc system: 

N r, 1 

( 3 - 6 ) W = £ -mig + Uixi) +'^2<> ij (\x i -x j \), 

t— 1 ^ J i<j 

where £, and x, arc the phase space of the i-the particle, m, is the particle mass, U(x,) is an external field 
and 4>tj(\Xi — Xj\) is a mean-field two-body interaction potential, the partition function is 

(3-7) Z = J dxd£exp(-H/ki t T) , 

where (x, £) represents the entire phase space of the /V-partic e system. Consequently, the free energy is 
given by 

(3.8) V = -kuTlnZ . 

Thus, information is neither gained nor lost whether the problem is formulated in terms of H or The 
advantage of using the free energy is that ^ is a global variable of state and therefore it is independent of 
coordinates. 
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4. Conclusion. In summary, we have carried out a systematic derivation of the lattice Boltzmann 
equation describing multi-phase flow from the Enskog equation a physically correct starting point for 
non-ideal gases. The model derived here is free of the defects of the existing models. The approach is 
rigorous and systematic. Not only the equation of state for non-ideal gases is obtained, but also the required 
thermodynamic consistency is achieved. Also, the procedure illustrated here is general and can be easily 
extended to other lattice Boltzmann models for complex fluids, e.g., binary mixtures or multi-component 
fluids. 
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